Physical observations of the region indicate that the NW Atlantic has experienced a regime shift that has accelerated warming in the area. Temperature has a profound impact on the growth, behavior, production, and mortality in ectotherms. For these reasons species are expected to follow their thermal preferences, otherwise risking consequences in the forms of metabolic costs and/or other forms of stress. Using trawl survey data we monitored we tracked the individual responses of 17 species during the ten years prior, and ten years following the regime shift in 2010. Species response to warming was studied in terms of their behavioral shifts via shifting geographic center of biomass,and their biologic response via their size-at-age relationships through time. The majority of species in the southern half of the study area shifted their center of biomass Northward, with no northern species shifting anywhere south of Cape Cod. A concurrent pattern emerged in growth patterns, with the majority of species experiencing faster early-stage growth and smaller adult size. These patterns suggest that even among species that are able to follow thermal preferences, there are additional metabolic costs associated with the shifting environment. Changes in growth rates across such a wide range of species have implications for ecologists and fisheries managers seeking to understand climate impacts on fisheries.
Introduction
Literature on TSR and Growth Expectations
Introduce Study Area and Existing Research
Recent changes in Gulf stream positioning have altered the relative influences of both the Gulf stream and Labrador current on temperature and salinity regimes for the region. Understanding that this new temperature and salinity regime should alter both the energetics of the region, as well as food-web dynamics. We seek to identify whether these changes have had a measurable impact on the growth of several groundfish species.
Despite many of these species living at or near the bottom, research suggests that the forces driving these thermal regimes should have far-reaching impacts biological impacts. These impacts may directly change the near-bottom environment through mixing patterns, and/or indirectly through food-web impacts.
Through the lens of an environmental regime change we will assess changes to body size, size-at-age, and growth characteristics derived using the Von-Bertalanffy growth equation.
Methods
Sea Surface Temperature
Global Sea surface temperature data was obtained via NOAA’s optimally interpolated SST analysis (OISSTv2), providing daily temperature values at a 0.25° latitude x 0.25° longitude resolution (Reynolds et al. 2007). A daily climatology for every 0.25° pixel in the global data set was created using average daily temperatures spanning the period of 1982-2011. Daily anomalies were then computed as the difference between observed temperatures and the daily climatological average. OISSTv2 data used in these analyses were provided by the NOAA PSL, Boulder, Colorado, USA from their website at https://psl.noaa.gov.
Temperature data was regionally averaged to match the sampling area for the fisheries independent survey program providing the age-at-length data. SST anomalies were averaged over the entire sampling region, consisting of continental shelf habitats from Cape Hatteras to Nova Scotia, to produce a daily time series. This time series was then processed into an annual timeseries of surface temperatures and anomalies. All area-averaging was done with area-weighting to account for differences in the relative areas of the latitude/longitude grid of the OISSTv2 data.
We computed the linear trend (\(°C/Decade\)) in the daily Gulf of Maine SST time series using standard linear regression. We considered the two time periods: all years before, and all years following the regime shift. Warming rates and long-term averages were estimated independently for each period.
Species Data
Fishery Independent data on groundfish size-at-age was collected as part of the NEFSC’s northeast trawl survey. This survey is conducted each year in the Spring and in the fall, with sample locations determined following a stratified-random survey design with effort allocated in proportion to stratum area. Add an Explanation of the Strata we Exclude. Correction factors were applied for changes in vessels, gear, and doors when appropriate. Prior to sampling, a CTD cast is performed to collect environmental conditions for the water column at/near the start of the trawl sample providing information on bottom habitat conditions like bottom temperatures.
Distribution Shift Analyses
Analyses focusing on the shifts in distribution among species were limited to data from the Spring survey season. Previous work has shown no significant changes in timing of sampling, or the mean annual latitude and longitude of the Spring survey sampling (Nye et al. 2009). Our analyses followed the continued movements of 30 species consistent with Nye et al. 2009. These species were originally selected for being consistently sampled across all years, and as representatives of a wide range of taxonomic groups.
Code
# Access the tidy survdat Data with {targets}withr::with_dir(rprojroot::find_root('_targets.R'), tar_load(survdat_clean))# Filter to Just spring?dist_dat <- survdat_clean %>%filter(season =="Spring")# Just the Species, or as Stocks?nye_species <-c("spiny dogfish","little skate","thorny skate","winter skate","alewife","american shad","atlantic herring","atlantic cod","haddock","pollock","silver hake","red hake","spotted hake","white hake","cusk","goosefish","american plaice","atlantic halibut","yellowtail flounder","winter flounder","fourspot flounder","summer flounder","windowpane","atlantic mackerel","ocean pout","atlantic wolffish","blackbelly rosefish","longhorn sculpin","sea raven","acadian redfish")# Filter to those speciesdist_dat <- dist_dat %>%filter(comname %in% nye_species)# Get the total biomass at each station# Totals up the abundance and biomass across sex*station_totals <- dist_dat %>%group_by(est_year, survey_area, stratum, tow, est_towdate, avgdepth, bottemp, decdeg_beglat, decdeg_beglon, comname) %>%summarise(biomass_kg =sum(biomass_kg, na.rm = T),.groups ="drop")
Calculating Spatial Metrics
Movements through time were characterized for each species using the following metrics: center of biomass, area-occupied, mean depth of occurrence, and the mean temperature of occurrence. Each metric was weighted by the biomass of a species sampled at each station such that.
\[X_j = \frac{ \sum_{}^{} w_iX_{ij} }{ \sum~w_i }\] Where \(X\) is the spatial metric, \(j\) is the year, & \(w\) is the biomass (in kg) for each station \(i\).
Code
# Function to run the group weighted means on different groups using ...# Does all variables: depth, temp, lat, lon & their variance (sd)group_spat_metrics <-function(station_totals, ...){ station_totals %>%group_by(comname, ...) %>%summarise(total_biomass =sum(biomass_kg),avg_biomass =mean(biomass_kg),biomass_sd =sd(biomass_kg),avg_depth =weightedMean(avgdepth, w = biomass_kg, na.rm = T),avg_temp =weightedMean(bottemp, w = biomass_kg, na.rm = T),avg_lat =weightedMean(decdeg_beglat, w = biomass_kg, na.rm = T),avg_lon =weightedMean(decdeg_beglon, w = biomass_kg, na.rm = T),depth_sd =weightedSd(avgdepth, w = biomass_kg, na.rm = T),temp_sd =weightedSd(bottemp, w = biomass_kg, na.rm = T),lat_sd =weightedSd(decdeg_beglat, w = biomass_kg, na.rm = T),lon_sd =weightedSd(decdeg_beglon, w = biomass_kg, na.rm = T),.groups ="drop")}# Function to calculate weighted mean & median across biomass variables# returns a common format that can join easily to global benchmarksbiomass_weighted_var <-function(station_totals, var_col, ...){ var_lab <-deparse(substitute(var_col)) mean_df <- station_totals %>%group_by(comname, ...) %>%summarise(biom_median =weightedMedian({{var_col}}, biomass_kg, na.rm = T),biom_mu =weightedMean({{var_col}}, biomass_kg, na.rm = T),biom_sd =weightedSd({{var_col}}, biomass_kg, na.rm = T), .groups ="drop") %>%mutate(var = var_lab)return(mean_df)}# 1. Yearly Spatial Metricsyear_avgs <-group_spat_metrics(station_totals, est_year) # 2. Average Metrics Across Temperature Regimesregime_avgs <- station_totals %>%filter(est_year >=2000) %>%mutate(temp_regime =ifelse(est_year <2010, "Early Regime", "Warm Regime")) %>%group_spat_metrics(., temp_regime) # Variance across all years within each speciesglobal_benchmarks <-group_spat_metrics(station_totals)
Analysis of Spatial Metrics
Size at Age Analyses
The available age data is a subset of the overall catch data from the NMFS trawl survey. This biological data subset contains additional information on individual fishes such as otolith-derived aging that require additional workup to ascertain. Not all fish species have age available for the full time period, with special attention being given to aid in the management of commercially valuable species. The data used for size-at-age analyses comes from both the Spring and Fall surveys.
Code
# Access Biological Data with {targets}withr::with_dir(rprojroot::find_root('_targets.R'), tar_load(survdat_biological))# 1. Prepare Bio Data for Age analysis for regimes:# Drop NA age data, restrict to summer and fallnmfs_bio <- survdat_biological %>%filter(is.na(age) ==FALSE, season %in%c("Spring", "Fall"))%>%as.data.frame()#Add regime and decade labelsnmfs_bio <- nmfs_bio %>%mutate(yearclass = est_year - (age-1),regime =ifelse(est_year <2010, "Early Regime", "Warm Regime"),decade =floor_decade(est_year)) rm(survdat_biological)#### Growth Data Prep #### # Filter to just the data for the two regimesregime_dat <-filter(nmfs_bio, est_year >=2000, est_year <=2019)# Rank species by how many measurements there arespecies_abunds <- regime_dat %>%count(comname) %>%arrange(desc(n)) # ordered by number measured# # Reorder alphabetically And take# vonbert_species <- sort(species_abunds$comname[1:17])# Drop the ones that don't workvonbert_species <- species_abunds %>%filter(comname %not in%c("goosefish", "fourspot flounder", "cusk", "ocean pout", "bluefish")) %>%arrange(comname) %>%pull(comname)
Size at age analyses were performed on 18 species. Species were omitted from analyses if age data was not available across both temperature regimes (bluefish & ocean pout) or if their physiology prevented accurate aging (elasmobranch species). The following 18 species had sufficient age data across both temperature regimes to be included in the analysis: acadian redfish, american plaice, atlantic cod, atlantic herring, atlantic mackerel, black sea bass, butterfish, haddock, pollock, red hake, scup, silver hake, summer flounder, white hake, windowpane, winter flounder, witch flounder and yellowtail flounder.
Code
# Name the list so it carries throughnames(vonbert_species) <- vonbert_species##### Reorganize species data into the list, and make group id for latervonbert_species_agedat <-map_dfr(vonbert_species, function(species_name){# Filter to the vonbert species regime_dat %>%filter(comname == species_name) %>%mutate(group_id =str_c(regime, comname, sep ="-")) }, .id ="common name") # Species we want panels for:key_species <-sort(c("american plaice","atlantic cod","atlantic herring","butterfish","haddock", "pollock","scup","silver hake","summer flounder","winter flounder","yellowtail flounder"))
Growth by Temperature Regime
To assess the impact of an elevated temperature regime, the size-at-age data was grouped into decadal periods relating to the shift in surface temperature regimes (2000-2009 & 2010-2019). Growth in fishes is commonly modeled using the “Von-Bertalanffy” growth function (VBGF) to capture how a fishes size (length or weight) changes with age. For each species the growth was modeled by fitting the VGBF to the length distribution-at-age.
\[L_t~=~L_{\infty}(1-e^{-K(t-t_0)})\] Where \(L_t\) is the length in cm of a species at age \(t\). \(L_\infty\) is the asymptotic maximum length, \(K\) is _____, and \(t_0\) is x intercept, or the hypothetical age at 0 cm length.
Code
# Set the function for nls() to solvevb <-vbFuns(param="Typical")
Code
# Function to Pull Vonbert Coefficients, # In development: boot strapped intervalsvonbert_coef <-function(length_age_dat, start_points, min_obs =20){#### Von Bert Fitting Using: {FSA} ##### Don't run on fewer than a minimum number of observations num_aged <-nrow(length_age_dat)# If number aged less than the minimum observations:# fill output with NA's, but preserve structure to build tableif(num_aged < min_obs){# dataframe structure for coefficients na_coef_df <-data.frame("Linf"=NA,"K"=NA,"t0"=NA)# dataframe structure for parameters na_param_df <-data.frame("parameter"=NA,"estimate"=NA,"std_error"=NA,"t_value"=NA,"pr_t"=NA,"converged"=NA)# What to return when it fails/aborts:return(list("model"=NA, "coef"= na_df, "params"= na_param_df)) } # Close conditional for sample size#### Fitting Model When Sample Size is Sufficient ##### Use nls to estimate VBGF parameters using starting points vbert_fit <-nls(length_cm ~vb(age, Linf, K, t0),data = length_age_dat,start = start_points)# Access parameters with coef() vbert_coef <-as.data.frame(t(coef(vbert_fit)))# Get Summary Info (For convergence and std. error) vbert_summ <-summary(vbert_fit)# Convergence Info vbert_conv <- vbert_summ$convInfo$isConv# Parameters w/ Standard Error vbert_params <-as.data.frame(vbert_summ$parameters) %>%rownames_to_column(var ="parameter") %>% janitor::clean_names()# add convergence information in vbert_params$converged <- vbert_conv#### Work in Progress / Debugging ##### This section fails, can't find namespace/object "length_age_dat"# # Use bootstrapping for confidence intervals:# vbert_boot <- Boot(vbert_fit) # Be patient! Be aware of some non-convergence# vbert_conf <- confint(vbert_boot)# output list out_l <-list("model"= vbert_fit,"coef"= vbert_coef,"params"= vbert_params#,#"conf" = vbert_conf )# Return the listreturn(out_l)}
Code
# Second Function for running that for a single Species, split by a factor column# Function does these things:# 1. grab data for one species# 2. Get starting points for VB using all data on that species# 3. Fit Vonbert to both Regimesprocess_group_vbert <-function(spec_name, split_col, min_obs){# Get the data using the species spec_dat <- dplyr::filter(vonbert_species_agedat, comname == spec_name)# Get starting points from all data test_starts <-vbStarts(length_cm ~ age, data = spec_dat)# Get coefficients for groups spec_dat %>%split(.[split_col]) %>%map(.f = vonbert_coef,start_points = test_starts, min_obs = min_obs) %>%# map_dfr(pluck, "coef", .id = split_col) %>% # Pull out just the coefficientsmap_dfr(pluck, "params", .id = split_col) %>%# Pull out the parameters and stderrmutate(comname = spec_name)}
Code
# Running Each Species * Regimespecies_coef <- vonbert_species %>%map_dfr(function(x){process_group_vbert(spec_name = x, split_col ="regime", min_obs =30) })
Growth Impacts on Productivity
Yield-per-recruit analysis similar to Baudron et al. 2014?
Results
Temperature Regime Shift
Code
# Pull the whole area's tempssurvey_temps <-filter(temp_regimes, survey_area =="Full Region")# Pull the regime averages outregime_avgs <- survey_temps %>%distinct(regime, avg_temp_c, avg_anom_c, temp_label) # Get the regime averages separately for in text valuesregime_1 <- regime_avgs %>%filter(regime =="1982-2009 Average")regime_2 <- regime_avgs %>%filter(regime =="2010-2021 Average")# Functions developed in 10_change_points.Rlibrary(bcp)# run bcp and flag changepoints using thresholdflag_bcp <-function(response_var, predictor_var =NULL, id_var, burnin =10000,mcmc =10000, return.mcmc = T, threshold_prob =0.3){# Use bcp to find support for change points bcp_x <-bcp(y = response_var, x = predictor_var, burnin = burnin, mcmc = mcmc, return.mcmc = return.mcmc)sink("/dev/null")# Convert summary to a dataframe bcp_sum <-as.data.frame(summary(bcp_x))sink()# Label the id variable bcp_sum$id <- id_var#Only pull records with support above desired threshold flag_dat <- bcp_sum[which(bcp_x$posterior.prob > threshold_prob), ] flag_dat <-as.data.frame(flag_dat)return(flag_dat)}# Set seedset.seed(3905)# Perform BCP regions temp breakpointsanom_bcp <-flag_bcp(response_var = survey_temps$sst_anom, id_var = survey_temps$yr, threshold_prob =0.3)# in-text valuesprob_text <-paste(round(anom_bcp$Probability *100, 2), collapse =", ")
SST anomaly time series of the region indicate a jump in annual temperature anomalies indicative of a regime shift centered around 2010. A bayesian change point analysis showed highest support for a regime shift among 3 the years: 2009-2011. Change point probabilities for these years were 33.66, 32.51, 32.49% respectively.
Code
# Function to plot bcp breaks, highlighting their supportplot_bcp_breaks <-function(bcp_data, bcp_breaks, region_label, y_label ="sizeSpectra Exponent (b)",y_col ="b"){# Set label elevation label_y <-mean(as.data.frame(bcp_data)[, y_col], na.rm = T) -1# plot with breakpointggplot(bcp_data, aes_string("yr", y_col)) +geom_line() +geom_point() +geom_vline(data = bcp_breaks, aes(xintercept = id), linetype =2, color ="gray50") + ggrepel::geom_label_repel(data = bcp_breaks, aes(x = id, y = label_y, label =str_c(round(Probability, 2), "%"))) +labs(x ="Year", y = y_label, title =str_c(region_label, " - Bayesian Change Points"))}# Plot the bcp breaks for the Sea surface temperature anomaliesplot_bcp_breaks(bcp_data = survey_temps, bcp_breaks = anom_bcp, y_label ="Temperature Anomalies", region_label ="Northeastern US Continental Shelf", y_col ="sst_anom")
Temperatures prior to the 2010 regime shift rose at a rate of 0.7 \(°C/Decade\). Regional temperatures then peaked in 2012 and have remained high while decreasing slowly at a rate of 0.02 \(°C/Decade\). Average temperatures following the regime shift (2011-2021), were on average 1.09°C warmer than the than the average over the previous 28 years (1982-2009). Temperatures within each regime were relatively stable, indicative of a jump in mean temperatures rather than a rate change between them. While temperatures in the current regime are falling, the rate is such that
Code
# Plot how each one exhibits a break in the temperature across the regiontemp_regimes %>%filter(survey_area =="Full Region") %>%ggplot( aes(yr, sst_anom)) +geom_col(aes(fill = regime), size =0.75, alpha =0.4) +scale_fill_gmri() +# Using geomtextpath for labelgeom_textpath(aes(x = yr, y = avg_anom_c, label = temp_label, linecolor = regime,colour = regime), size =4, vjust =-1, fontface ="bold", show.legend =FALSE) +# Line segments for the regime averages:geom_segment(aes(x = x, xend = xend,y = avg_anom_c, yend = avg_anom_c,color = regime),#key_glyph = "rect",size = .8,linetype =1) +scale_x_continuous(expand =expansion(add =c(0.25,0.25))) +scale_y_continuous(labels =number_format(suffix =" \u00b0C")) +scale_color_gmri() +guides(fill ="none") +theme(legend.title =element_blank(),legend.position ="bottom",legend.background =element_rect(fill ="transparent"),legend.key =element_rect(fill ="transparent", color ="transparent")) +labs(title ="Northwest-Atlantic Temperature Regime Shift",x ="", y ="Surface Temperature Anomaly",caption ="Anomalies calculated using 1982-2011 reference period.") +guides(label ="none",color =guide_legend(keyheight =unit(0.5, "cm"), ))
Distribution Shifts
Code
# Get the values for each of them in a standardized format# Standardize to global benchmarks# Latitudelat_yrly <-biomass_weighted_var(station_totals, decdeg_beglat, est_year) %>%left_join(select(global_benchmarks, c("comname", "avg_lat", "lat_sd")), by ="comname") %>%mutate(biom_z = (biom_mu - avg_lat) / lat_sd) # Longitudelon_yrly <-biomass_weighted_var(station_totals, decdeg_beglon, est_year) %>%left_join(select(global_benchmarks, c("comname", "avg_lon", "lon_sd")), by ="comname") %>%mutate(biom_z = (biom_mu - avg_lon) / lon_sd) # Bottom Temperaturet_yrly <-biomass_weighted_var(station_totals, bottemp, est_year) %>%left_join(select(global_benchmarks, c("comname", "avg_temp", "temp_sd")), by ="comname") %>%mutate(biom_z = (biom_mu - avg_temp) / temp_sd) # Depthd_yrly <-biomass_weighted_var(station_totals, avgdepth, est_year) %>%left_join(select(global_benchmarks, c("comname", "avg_depth", "depth_sd")), by ="comname") %>%mutate(biom_z = (biom_mu - avg_depth) / depth_sd) # Put them together in one table for heatmapstandardized_changes <-list("Latitude"= lat_yrly,"Longitude"= lon_yrly,"Bottom Temperature"= t_yrly,"Depth"= d_yrly) %>%map_dfr(function(x){select(x, comname, est_year, var, biom_z)}, .id ="var_neat")# Organize them by their global average latitudeglobal_benchmarks <- global_benchmarks %>%mutate(comname =fct_reorder(.f = comname, .x = avg_lat, .fun = median))# Relevel the standardized changesstandardized_changes <- standardized_changes %>%mutate(comname =factor(comname, levels =levels(global_benchmarks$comname)))# Heatmapstandardized_changes %>%ggplot(aes(est_year, comname, fill = biom_z)) +geom_tile() +scale_x_continuous(expand =expansion(add =c(0, 0))) +facet_wrap(~var_neat, ncol =2, scales ="fixed") +scale_fill_distiller(palette ="RdBu",oob = scales::oob_squish,limits =c(-1, 1), breaks =c(-1, -.5, 0, .5, 1),labels =c("<-1", "-.5", "0", ".5", ">1")) +labs(x ="Year", y ="",title ="Center of Biomass Movements Through Time",fill ="Center of Biomass Metric Change (z)",caption ="Species ordered based on center of biomass latitude across all years.") +theme(legend.position ="bottom",axis.text.y =element_text(size =7)) +guides(fill =guide_colorbar(title.position ="top", title.hjust =0.5, barwidth =unit(6, "cm")))
# Standardized Plotslat_yrly %>%ggplot(aes(est_year, biom_z, group = comname)) +# geom_path(aes(group = comname), linetype = 1, color = "gray50", size = 0.5, alpha = 0.2) +geom_smooth(method ="lm", formula = y ~ x, alpha =0.8, se =FALSE, color ="gray50") +geom_jitter(width =0.1, size =2, alpha =0.25) +geom_hline(aes(yintercept =0, color ="Overall Average"), size =1, linetype =1) +scale_color_gmri(reverse = T) +labs(title ="Shift in Center of Biomass Latitudes",subtitle ="Yearly averages weighted by biomass, scaled against all-year average.",y ="Center of Biomass Latitude (z)", x ="", color ="Time Period", fill ="Time Period") +theme(legend.position ="bottom")
Code
# Longitudelon_yrly %>%ggplot(aes(est_year, biom_z, group = comname)) +# geom_path(aes(group = comname), linetype = 1, color = "gray50", size = 0.5, alpha = 0.2) +geom_smooth(method ="lm", formula = y ~ x, alpha =0.8, se =FALSE, color ="gray50") +geom_jitter(width =0.1, size =2, alpha =0.25) +geom_hline(aes(yintercept =0, color ="Overall Average"), size =1, linetype =1) +scale_color_gmri(reverse = T) +#coord_flip() +labs(title ="Shift in Center of Biomass Longtitudes",subtitle ="Yearly averages weighted by biomass, scaled against all-year average.",y ="Center of Biomass Longitude (z)", x ="", color ="Time Period", fill ="Time Period") +theme(strip.text.y =element_text(size =7, angle =0, hjust =0),legend.position ="bottom")
Code
# Temperaturet_yrly %>%ggplot(aes(est_year, biom_z, group = comname)) +# geom_path(aes(group = comname), linetype = 1, color = "gray50", size = 0.5, alpha = 0.2) +geom_smooth(method ="lm", formula = y ~ x, alpha =0.8, se =FALSE, color ="gray50") +geom_jitter(width =0.1, size =2, alpha =0.25) +geom_hline(aes(yintercept =0, color ="Overall Average"), size =1, linetype =1) +scale_color_gmri(reverse = T) +#coord_flip() +labs(title ="Shift in Center of Biomass Temperature",subtitle ="Yearly averages weighted by biomass, scaled against all-year average.",y ="Center of Biomass Bottom Temperature (z)", x ="", color ="Time Period", fill ="Time Period") +theme(legend.position ="bottom")
Code
# Depthd_yrly %>%ggplot(aes(est_year, biom_z, group = comname)) +# geom_path(aes(group = comname), linetype = 1, color = "gray50", size = 0.5, alpha = 0.2) +geom_smooth(method ="lm", formula = y ~ x, alpha =0.8, se =FALSE, color ="gray50") +geom_jitter(width =0.1, size =2, alpha =0.25) +geom_hline(aes(yintercept =0, color ="Overall Average"), size =1, linetype =1) +scale_color_gmri(reverse = T) +#coord_flip() +labs(title ="Shift in Center of Biomass Depth",subtitle ="Yearly averages weighted by biomass, scaled against all-year average.",y ="Center of Biomass Depth (z)", x ="", color ="Time Period", fill ="Time Period") +theme(strip.text.y =element_text(size =7, angle =0, hjust =0),legend.position ="bottom")
Changes to Growth
Code
# Take the first part from the previous plot:# takes each species# for age 0 - max_age, solve what the size would be based on vb fitregime_curve_prep <-function(spec_name){# The raw size/age data: t_dat <- vonbert_species_agedat %>%filter(comname == spec_name) # The group coefficients, augmented with x values across a range age_max <-max(t_dat$age) x_range <-seq(from =0, to = age_max, by = .1) t_coef <- species_coef %>%distinct(comname, regime, parameter, estimate) %>%pivot_wider(values_from ="estimate", names_from ="parameter") %>%filter(comname == spec_name) %>%right_join(# Make a dataframe that covers the size range for the speciesy =data.frame(regime =rep(unique(t_dat$regime), length(x_range)),age =rep(x_range, length(unique(t_dat$regime)))),by ="regime") %>%mutate(# Calculate the length at agepred_length = Linf * (1-exp(-1* K * (age - t0))) )# Pivot wider?# Take the size at age from each regime and get the difference t_coef <- t_coef %>%select(-c(Linf, K, t0)) %>%pivot_wider(names_from ="regime", values_from ="pred_length") %>%mutate(growth_shift =`Warm Regime`-`Early Regime`,growth_frac = (growth_shift /`Early Regime`),size_change =ifelse(growth_shift >0, "Increase", "Decrease") )return(t_coef)}# Do it for all the speciesregime_comparison <- vonbert_species %>%map_dfr(.f = regime_curve_prep, .id ="comname")
Changes in the size-at-age growth relationship varied across species, with the majority (14/18) of the species exhibiting smaller L-infinity during the warmer temperature regime.
Code
# Try geom_waffle to just count each species as an individual and show that they are growing to smaller sizes in new regime:library(ggwaffle)age_set <-4ggplot(data =waffle_iron(data =filter(regime_comparison, age == age_set), aes_d(group = size_change), rows =2), aes(x, y, fill =fct_rev(group))) +geom_waffle() +theme_waffle() +scale_fill_gmri() +theme(axis.text =element_blank(),axis.text.x =element_blank(),axis.text.y =element_blank(),legend.title =element_text(face ="bold", size =12),legend.position ="bottom",panel.grid.major.y =element_blank()) +labs(x ="", y ="", fill =str_c("Warm Regime Growth Impacts on ", age_set," Body Length:"))
The species that showed larger asymptotic lengths include the three hake species: red hake, silver hake, & white hake as well as Atlantic herring. These species are all among the majority group of species that is found primarily in the Northern part of the survey area, in and around the Gulf of Maine.
Code
# Reshape things so you can show confidence intervals hereplot_dat <- species_coef %>%filter(parameter !="t0") %>%mutate(llim = estimate -2* std_error,rlim = estimate +2* std_error,residence =ifelse( comname %in%c("atlantic mackerel", "scup", "black sea bass", "butterfish"),"Southern\nResident", "Gulf of Maine\nResident"),regime =fct_rev(regime),comname =fct_rev(comname))# Plot itggplot(plot_dat, aes(x = estimate, y = comname, group = regime, color = regime)) +# need to nudge groups separately# Early Regimegeom_segment(data =filter(plot_dat, regime =="Early Regime"), aes(x = llim, xend = rlim, yend = comname), position =position_nudge(x =0, y =0.25)) +geom_segment(data=filter(plot_dat, regime =="Warm Regime"), aes(x = llim, xend = rlim, yend = comname), position =position_nudge(x =0, y =-0.25)) +geom_point(aes(color = regime), position =position_dodge(width =1), size =2) +#facet_grid(. ~ parameter, scales = "free", space = "free_y") +facet_grid(residence ~ parameter, scales ="free", space ="free_y") +scale_color_gmri(reverse = F) +labs(x ="Coefficient Estimate", y ="", fill ="Temperature Regime:") +theme(legend.position ="bottom", panel.background =element_rect(color ="gray80", fill ="transparent", size =1),panel.spacing =unit(2, "lines"),panel.border =element_rect(color ="gray80", fill ="transparent", size =1))
Code
# Plotting Everythingspecies_coef %>%distinct(comname, regime, parameter, estimate) %>%pivot_wider(values_from =c("estimate"), names_from ="parameter") %>%pivot_longer(names_to ="vbcoef", values_to ="val", cols =c("K", "Linf", "t0")) %>%filter(vbcoef !="t0") %>%mutate(residence =ifelse(comname %in%c("atlantic mackerel", "scup", "black sea bass", "butterfish"), "Southern\nResident", "Gulf of Maine\nResident")) %>%ggplot(aes(x = val, y =fct_rev(comname), fill =fct_rev(regime))) +geom_col(position ="dodge") +# facet_grid(. ~ vbcoef, scales = "free", space = "free_y") +facet_grid(residence ~ vbcoef, scales ="free", space ="free_y") +scale_fill_gmri(reverse = T) +labs(x ="Coefficient Estimate", y ="", fill ="Temperature Regime:") +theme(legend.position ="bottom", panel.background =element_rect(color ="gray80", fill ="transparent", size =1),panel.spacing =unit(2, "lines"),panel.border =element_rect(color ="gray80", fill ="transparent", size =1))
# Make the plot for each speciesvb_fits_plot =function(spec_name){# The raw size/age data: t_dat <- vonbert_species_agedat %>%filter(comname == spec_name) # The group coefficients, augmented with x values across a range age_max <-max(t_dat$age) x_range <-seq(from =0, to = age_max, by = .1) t_coef <- species_coef %>%distinct(comname, regime, parameter, estimate) %>%pivot_wider(values_from ="estimate", names_from ="parameter") %>%filter(comname == spec_name) %>%right_join(y =data.frame(regime =rep(unique(t_dat$regime), length(x_range)),age =rep(x_range, length(unique(t_dat$regime)))),by ="regime") %>%mutate(# Calculate the length at agepred_length = Linf * (1-exp(-1* K * (age - t0))) )# Make plot?ggplot(t_dat, aes(x = age, y = length_cm)) +geom_point(alpha =0.25) +#facet_wrap(~regime, ncol = 1) +geom_line(data = t_coef,aes(x = age, y = pred_length, group = regime, color = regime), size =1) +scale_x_continuous(breaks =seq(0, age_max, by =2),limits =c(0, age_max),expand =expansion(mult =c(0.1, .2))) +scale_color_gmri() +geom_dl(data = t_coef, aes(x = age, y = pred_length, label = regime, color = regime), method=list("last.qp", cex =0.65, hjust =-.15),size =12) +labs(x ="Age",y ="Length (cm)",title = spec_name,color ="Temperature Regime:",subtitle ="Length at Age Regimes") }
Code
# Loop through the key species:walk(key_species, plot_panelset, plot_fun = vb_fits_plot)
american plaice
atlantic cod
atlantic herring
butterfish
haddock
pollock
scup
silver hake
summer flounder
winter flounder
yellowtail flounder
Discussion
Bonus: Data Explorations
Age Data Availability
The availability of age data varied by species and through time, with older ages becoming less frequently sampled in recent years across many species.
Code
# Plotting age distribution as ggridgesplot_age_ridgeplot <-function(species){ vonbert_species_agedat %>%filter(comname == species, age <=15) %>%mutate(yr =factor(est_year),yr =fct_rev(yr)) %>%ggplot(aes(x = age, y = yr)) +geom_density_ridges(aes(fill = regime), color ="white", show.legend =FALSE) +scale_fill_gmri() +facet_wrap(~comname) +scale_x_continuous(limits =c(0,NA), expand =expansion(add =c(0, 0))) +labs(x ="Age", y ="Year") +guides(size =guide_legend(title.position ="top", title.hjust =0.5))}
Code
# Loop through the key species:walk(key_species, plot_panelset, plot_fun = plot_age_ridgeplot)
american plaice
atlantic cod
atlantic herring
butterfish
haddock
pollock
scup
silver hake
summer flounder
winter flounder
yellowtail flounder
Strong age classes were particularly prominent in the following species, showing lasting numbers in subsequent years.
Species
Strong Year Classes
Atlantic Cod
2003, 2013
Haddock
2003, 2010, 2013
Atlantic Herring
2008, 2011
Pollock
2000, 2012
Growth Increment Changes
Annual Growth increments for each species were computed to track the change in growth of specific age-groups through time.
Code
# Age-Size Increments Differencesvb_incs <- vonbert_species_agedat %>%group_by(`common name`, est_year, age, regime) %>%summarise(avg_len =mean(length_cm, na.rm = T),avg_wt =mean(indwt, na.rm = T),.groups ="drop") %>%arrange(`common name`, age) %>%group_split(`common name`, est_year) %>%map_dfr(function(x){ x %>%mutate(# What age is being subtracted fromlead_age =lead(age, n =1),# String built from the value to confirmdiff_name =str_c(age, " - ", lead_age),# Length/Weight at age + 1len_age_p1 =lead(avg_len, n =1),wt_age_p1 =lead(avg_wt, n =1),# Length/Weight Growth Incrementslen_inc = len_age_p1 - avg_len,wt_inc = wt_age_p1 - avg_wt,# Percent change in bodymass for going from age to age+1len_perc_change = (len_inc / avg_len) *100,wt_perc_change = (wt_inc / avg_wt) *100) %>%filter((lead_age - age ==1)) })
Age specific growth (annual growth increments) help detail where in a species life cycle it may be experiencing favorable or less-favorable conditions for growth. These can be a result of differences in habitat use or prey resources at different ages, and/or a change in exposure to temperature or other stressors that inflict a metabolic cost.
Code
# Put together a timeseries of:# how big were the age 0's,# how big were the age 1's# Plotting Age Incrementsplot_age_increments <-function( spec_choice, var_choice ="length", gam_knots =15, facet_prefix ="Age Increment: ",facet_suffix =""){# Code to switch variables and their labels: var_sym <-switch(var_choice,"length"=sym("len_inc"),"weight"=sym("wt_inc")) var_label <-switch(var_choice,"length"="Annual Growth (cm)","weight"="Annual Growth (kg)")# Add pre/suffix to facet labels appender <-function(string, prefix = facet_prefix, suffix = facet_suffix) {paste0(prefix, string)}# Make the plot vb_incs %>%filter(`common name`== spec_choice, lead_age <=5) %>%drop_na({{var_sym}}) %>%ggplot(aes(est_year, {{var_sym}})) +# geom_col(fill = gmri_cols("gmri blue")) +geom_point(color =gmri_cols("gmri blue")) +geom_line(color =gmri_cols("gmri blue"), group =1, alpha =0.3, size =0.8) +geom_smooth(method ="gam",formula = y ~s(x, bs ="cs", k = gam_knots),se =FALSE,size =1,color =gmri_cols("orange")) +facet_wrap(~diff_name, ncol =1, scales ="free",labeller =as_labeller(appender)) +theme(legend.position ="bottom") +labs(x ="",y = var_label,color ="Age Increment",title =str_to_title(spec_choice),subtitle =str_c("Changes in annual ", var_choice, " growth increments."))}
Code
# Loop through the key species:walk(key_species, plot_panelset, plot_fun = plot_age_increments, var_choice ="weight", gam_knots =5)
american plaice
atlantic cod
atlantic herring
butterfish
haddock
pollock
scup
silver hake
summer flounder
winter flounder
yellowtail flounder
Length at Age Changes
Code
# Distribution of Length/Weight Across Regimes plot_regime_sizes <-function(spec, var ="length", age_lim =5, facet_prefix ="Age - ", facet_suffix =""){# switches for length or width var_col <-switch (var,"length"=sym("length_cm"),"weight"=sym("indwt")) var_lab <-switch (var,"length"="Individual Length (cm)","weight"="Individual Weight (kg)")# Add pre/suffix to facet labels appender <-function(string, prefix = facet_prefix, suffix = facet_suffix) {paste0(prefix, string)}# Select Ages, Prep data size_changes <- regime_dat %>%filter( comname == spec, age <=6) %>%## Calculate group averagesgroup_by(age, regime) %>%mutate(median_size =median(!!{{var_col}})) %>%ungroup() # Make median dataframe median_df <- size_changes %>%distinct(age, regime, median_size) %>%pivot_wider(names_from ="regime", values_from ="median_size")# Make Plot size_changes %>%ggplot(aes(x =!!{{var_col}})) +## distributionstat_halfeye(aes(color = regime, fill =after_scale(color)),slab_alpha = .5, .width =0, trim = T, shape =21, point_colour ="gray25",stroke =1.6,scale = .86) +# median linegeom_segment(data = median_df,aes(x =`Early Regime`,xend =`Warm Regime`,y =0,yend =0),size = .8,color ="grey25", # , position = position_nudge(y = -.15) ) +## median points - not workingstat_halfeye(aes(color = regime), .width =c(0), slab_fill =NA) +facet_wrap(~age, ncol =1, labeller =as_labeller(appender), scales ="free") +scale_color_gmri() +scale_fill_gmri() +scale_y_continuous(labels =percent_format()) +theme(legend.position ="bottom") +labs(x = var_lab,y ="Proportion",title = spec,fill ="Temperature Regime",color ="Temperature Regime",shape ="Median" )}# Tester: # plot_regime_sizes(key_species[[1]], var = "length")
Code
# Loop through the key species:walk(key_species, plot_panelset, plot_fun = plot_regime_sizes, var ="length")
american plaice
atlantic cod
atlantic herring
butterfish
haddock
pollock
scup
silver hake
summer flounder
winter flounder
yellowtail flounder
Weight at Age Changes
Code
# Loop through the key species:walk(key_species, plot_panelset, plot_fun = plot_regime_sizes, var ="weight")
american plaice
atlantic cod
atlantic herring
butterfish
haddock
pollock
scup
silver hake
summer flounder
winter flounder
yellowtail flounder
Shifts in Growth
Code
# Plot their growth curvesregime_comparison %>%#filter(comname %in% key_species) %>% ggplot(aes(x = age, y = growth_frac)) +geom_hline(yintercept =0, size =0.5, lty =2, color ="black") +geom_line(aes(group = comname), show.legend =FALSE, size =1, color ="orange")+scale_y_continuous(labels =percent_format()) +facet_wrap(~comname, scales ="free_x", ncol =4) +labs(x ="Age", y ="Shift in Length Moving to Warm Regime", title ="Thermal Regime Impacts on Growth",subtitle ="Changes in Size at Age (Warm Regime - Historic Regime)")
Code
# Kathy requested the curve data.# Going to save both the coefficients and the fitted data for all specieslibrary(here)write_csv(regime_comparison, here("data/growth_results/growth_regime_fits.csv"))write_csv(species_coef, here("data/growth_results/growth_regime_coefs.csv"))
Erdman, Chandra, and John Emerson. 2007. “Bcp: An r Package for Performing a Bayesian Analysis of Change Point Problems.”Journal of Statistical Software 23 (December). https://doi.org/10.18637/jss.v023.i03.
Nye, Ja, Js Link, Ja Hare, and Wj Overholtz. 2009. “Changing Spatial Distribution of Fish Stocks in Relation to Climate and Population Size on the Northeast United States Continental Shelf.”Marine Ecology Progress Series 393 (October): 111–29. https://doi.org/10.3354/meps08220.
Reynolds, Richard W., Thomas M. Smith, Chunying Liu, Dudley B. Chelton, Kenneth S. Casey, and Michael G. Schlax. 2007. “Daily High-Resolution-Blended Analyses for Sea Surface Temperature.”Journal of Climate 20 (22): 5473–96. https://doi.org/10.1175/2007JCLI1824.1.
Source Code
---title: "Growth Regimes in the NW Atlantic Following Temperature Regime Shift"author: name: "Adam Kemberling" url: https://github.com/adamkemberling affiliation: Gulf of Maine Research Institutedescription: | Differential Growth Rates of Northeast US Fishes Following Temperature Regime Shiftdate: "`r Sys.Date()`"format: html: self-contained: true code-fold: true code-tools: true df-print: kable toc: true toc-depth: 2editor: sourceexecute: echo: true warning: false message: false fig.height: 6 fig.width: 6 fig.align: "center" comment: ""bibliography: references.bib---```{r packages and functions}#| include = FALSE,#| echo = FALSE#### Packages ####library(FSAdata)library(FSA)library(car)library(matrixStats)library(scales)library(patchwork)library(directlabels)library(xaringanExtra)library(ggridges)library(ggdist)library(ggforce)library(geomtextpath)library(colorspace)library(gt)library(tidyverse)library(gmRi)library(targets)# Set themetheme_set(theme_gmri())# Acitvate panelsetxaringanExtra::use_panelset()# Function to generate panels using a key word and a plot functionplot_panelset <-function(spec, plot_fun, ...) {# Open panelcat("::: {.panel}\n")# Create header for panel namecat("##", spec, "{.panel-name}\n")# # Text as a plot description# cat("Plot of", spec, "\n")# Make plot p <-plot_fun(spec, ...)print(p)# Add linebreakcat("\n")# Close panelcat(":::\n")}````r gmRi::use_gmri_style_rmd(css_file = "gmri_rmarkdown.css")`# AbstractPhysical observations of the region indicate that the NW Atlantic has experienced a regime shift that has accelerated warming in the area. Temperature has a profound impact on the growth, behavior, production, and mortality in ectotherms. For these reasons species are expected to follow their thermal preferences, otherwise risking consequences in the forms of metabolic costs and/or other forms of stress. Using trawl survey data we monitored we tracked the individual responses of 17 species during the ten years prior, and ten years following the regime shift in 2010. Species response to warming was studied in terms of their behavioral shifts via shifting geographic center of biomass,and their biologic response via their size-at-age relationships through time. The majority of species in the southern half of the study area shifted their center of biomass Northward, with no northern species shifting anywhere south of Cape Cod. A concurrent pattern emerged in growth patterns, with the majority of species experiencing faster early-stage growth and smaller adult size. These patterns suggest that even among species that are able to follow thermal preferences, there are additional metabolic costs associated with the shifting environment. Changes in growth rates across such a wide range of species have implications for ecologists and fisheries managers seeking to understand climate impacts on fisheries.# Introduction### Literature on TSR and Growth Expectations### Introduce Study Area and Existing ResearchRecent changes in Gulf stream positioning have altered the relative influences of both the Gulf stream and Labrador current on temperature and salinity regimes for the region. Understanding that this new temperature and salinity regime should alter both the energetics of the region, as well as food-web dynamics. We seek to identify whether these changes have had a measurable impact on the growth of several groundfish species.Despite many of these species living at or near the bottom, research suggests that the forces driving these thermal regimes should have far-reaching impacts biological impacts. These impacts may directly change the near-bottom environment through mixing patterns, and/or indirectly through food-web impacts.Through the lens of an environmental regime change we will assess changes to body size, size-at-age, and growth characteristics derived using the Von-Bertalanffy growth equation.# Methods## Sea Surface TemperatureGlobal Sea surface temperature data was obtained via NOAA's optimally interpolated SST analysis (OISSTv2), providing daily temperature values at a 0.25° latitude x 0.25° longitude resolution [@reynolds2007]. A daily climatology for every 0.25° pixel in the global data set was created using average daily temperatures spanning the period of 1982-2011. Daily anomalies were then computed as the difference between observed temperatures and the daily climatological average. OISSTv2 data used in these analyses were provided by the NOAA PSL, Boulder, Colorado, USA from their website at https://psl.noaa.gov.```{r}# Load the regional temperatures from {targets}withr::with_dir(rprojroot::find_root('_targets.R'), tar_load(regional_oisst))# Get regional averagestemp_regimes <- regional_oisst %>%filter(yr >1981) %>%mutate(regime =ifelse(yr <2010, "1982-2009 Average", "2010-2021 Average"),survey_area =ifelse(survey_area =="all", "Full Region", survey_area),survey_area =factor(survey_area, levels =c("Full Region", "GoM", "GB", "SNE", "MAB"))) %>%group_by(survey_area, regime) %>%mutate(avg_temp_c =round(mean(sst), 2),temp_label =str_c(avg_temp_c, "\u00b0C"),avg_anom_c =mean(sst_anom),#avg_anom_c = ifelse(regime == "1982-2009 Regime", 0, avg_anom_c),x =ifelse(regime =="1982-2009 Average", 1981.5, 2009.5),xend =ifelse(regime =="1982-2009 Average", 2009.5, 2021.5), ) %>%ungroup()```### Temperature Study RegionTemperature data was regionally averaged to match the sampling area for the fisheries independent survey program providing the age-at-length data. SST anomalies were averaged over the entire sampling region, consisting of continental shelf habitats from Cape Hatteras to Nova Scotia, to produce a daily time series. This time series was then processed into an annual timeseries of surface temperatures and anomalies. All area-averaging was done with area-weighting to account for differences in the relative areas of the latitude/longitude grid of the OISSTv2 data.### Temperature Regime ShiftAnnual temperature anomalies were checked for a regime shift using the bayesian change points package [@erdman2007] following the approach of [@Barry2012] implemented in the R statistical programming language.### Warming TrendsWe computed the linear trend ($°C/Decade$) in the daily Gulf of Maine SST time series using standard linear regression. We considered the two time periods: all years before, and all years following the regime shift. Warming rates and long-term averages were estimated independently for each period.## Species DataFishery Independent data on groundfish size-at-age was collected as part of the NEFSC's northeast trawl survey. This survey is conducted each year in the Spring and in the fall, with sample locations determined following a stratified-random survey design with effort allocated in proportion to stratum area. **Add an Explanation of the Strata we Exclude**. Correction factors were applied for changes in vessels, gear, and doors when appropriate. Prior to sampling, a CTD cast is performed to collect environmental conditions for the water column at/near the start of the trawl sample providing information on bottom habitat conditions like bottom temperatures.## Distribution Shift AnalysesAnalyses focusing on the shifts in distribution among species were limited to data from the Spring survey season. Previous work has shown no significant changes in timing of sampling, or the mean annual latitude and longitude of the Spring survey sampling [@nye2009]. Our analyses followed the continued movements of 30 species consistent with Nye et al. 2009. These species were originally selected for being consistently sampled across all years, and as representatives of a wide range of taxonomic groups.```{r}# Access the tidy survdat Data with {targets}withr::with_dir(rprojroot::find_root('_targets.R'), tar_load(survdat_clean))# Filter to Just spring?dist_dat <- survdat_clean %>%filter(season =="Spring")# Just the Species, or as Stocks?nye_species <-c("spiny dogfish","little skate","thorny skate","winter skate","alewife","american shad","atlantic herring","atlantic cod","haddock","pollock","silver hake","red hake","spotted hake","white hake","cusk","goosefish","american plaice","atlantic halibut","yellowtail flounder","winter flounder","fourspot flounder","summer flounder","windowpane","atlantic mackerel","ocean pout","atlantic wolffish","blackbelly rosefish","longhorn sculpin","sea raven","acadian redfish")# Filter to those speciesdist_dat <- dist_dat %>%filter(comname %in% nye_species)# Get the total biomass at each station# Totals up the abundance and biomass across sex*station_totals <- dist_dat %>%group_by(est_year, survey_area, stratum, tow, est_towdate, avgdepth, bottemp, decdeg_beglat, decdeg_beglon, comname) %>%summarise(biomass_kg =sum(biomass_kg, na.rm = T),.groups ="drop")```### Calculating Spatial MetricsMovements through time were characterized for each species using the following metrics: center of biomass, area-occupied, mean depth of occurrence, and the mean temperature of occurrence. Each metric was weighted by the biomass of a species sampled at each station such that.$$X_j = \frac{ \sum_{}^{} w_iX_{ij} }{ \sum~w_i }$$ Where $X$ is the spatial metric, $j$ is the year, & $w$ is the biomass (in kg) for each station $i$.```{r}# Function to run the group weighted means on different groups using ...# Does all variables: depth, temp, lat, lon & their variance (sd)group_spat_metrics <-function(station_totals, ...){ station_totals %>%group_by(comname, ...) %>%summarise(total_biomass =sum(biomass_kg),avg_biomass =mean(biomass_kg),biomass_sd =sd(biomass_kg),avg_depth =weightedMean(avgdepth, w = biomass_kg, na.rm = T),avg_temp =weightedMean(bottemp, w = biomass_kg, na.rm = T),avg_lat =weightedMean(decdeg_beglat, w = biomass_kg, na.rm = T),avg_lon =weightedMean(decdeg_beglon, w = biomass_kg, na.rm = T),depth_sd =weightedSd(avgdepth, w = biomass_kg, na.rm = T),temp_sd =weightedSd(bottemp, w = biomass_kg, na.rm = T),lat_sd =weightedSd(decdeg_beglat, w = biomass_kg, na.rm = T),lon_sd =weightedSd(decdeg_beglon, w = biomass_kg, na.rm = T),.groups ="drop")}# Function to calculate weighted mean & median across biomass variables# returns a common format that can join easily to global benchmarksbiomass_weighted_var <-function(station_totals, var_col, ...){ var_lab <-deparse(substitute(var_col)) mean_df <- station_totals %>%group_by(comname, ...) %>%summarise(biom_median =weightedMedian({{var_col}}, biomass_kg, na.rm = T),biom_mu =weightedMean({{var_col}}, biomass_kg, na.rm = T),biom_sd =weightedSd({{var_col}}, biomass_kg, na.rm = T), .groups ="drop") %>%mutate(var = var_lab)return(mean_df)}# 1. Yearly Spatial Metricsyear_avgs <-group_spat_metrics(station_totals, est_year) # 2. Average Metrics Across Temperature Regimesregime_avgs <- station_totals %>%filter(est_year >=2000) %>%mutate(temp_regime =ifelse(est_year <2010, "Early Regime", "Warm Regime")) %>%group_spat_metrics(., temp_regime) # Variance across all years within each speciesglobal_benchmarks <-group_spat_metrics(station_totals) ```### Analysis of Spatial Metrics## Size at Age AnalysesThe available age data is a subset of the overall catch data from the NMFS trawl survey. This biological data subset contains additional information on individual fishes such as otolith-derived aging that require additional workup to ascertain. Not all fish species have age available for the full time period, with special attention being given to aid in the management of commercially valuable species. The data used for size-at-age analyses comes from both the Spring and Fall surveys.```{r bio data prep}# Access Biological Data with {targets}withr::with_dir(rprojroot::find_root('_targets.R'), tar_load(survdat_biological))# 1. Prepare Bio Data for Age analysis for regimes:# Drop NA age data, restrict to summer and fallnmfs_bio <- survdat_biological %>%filter(is.na(age) ==FALSE, season %in%c("Spring", "Fall"))%>%as.data.frame()#Add regime and decade labelsnmfs_bio <- nmfs_bio %>%mutate(yearclass = est_year - (age-1),regime =ifelse(est_year <2010, "Early Regime", "Warm Regime"),decade =floor_decade(est_year)) rm(survdat_biological)#### Growth Data Prep #### # Filter to just the data for the two regimesregime_dat <-filter(nmfs_bio, est_year >=2000, est_year <=2019)# Rank species by how many measurements there arespecies_abunds <- regime_dat %>%count(comname) %>%arrange(desc(n)) # ordered by number measured# # Reorder alphabetically And take# vonbert_species <- sort(species_abunds$comname[1:17])# Drop the ones that don't workvonbert_species <- species_abunds %>%filter(comname %not in%c("goosefish", "fourspot flounder", "cusk", "ocean pout", "bluefish")) %>%arrange(comname) %>%pull(comname)```Size at age analyses were performed on `r length(vonbert_species)` species. Species were omitted from analyses if age data was not available across both temperature regimes (bluefish & ocean pout) or if their physiology prevented accurate aging (elasmobranch species). The following 18 species had sufficient age data across both temperature regimes to be included in the analysis: **`r pander::pander(vonbert_species)`**.```{r panel species}# Name the list so it carries throughnames(vonbert_species) <- vonbert_species##### Reorganize species data into the list, and make group id for latervonbert_species_agedat <-map_dfr(vonbert_species, function(species_name){# Filter to the vonbert species regime_dat %>%filter(comname == species_name) %>%mutate(group_id =str_c(regime, comname, sep ="-")) }, .id ="common name") # Species we want panels for:key_species <-sort(c("american plaice","atlantic cod","atlantic herring","butterfish","haddock", "pollock","scup","silver hake","summer flounder","winter flounder","yellowtail flounder"))```### Growth by Temperature RegimeTo assess the impact of an elevated temperature regime, the size-at-age data was grouped into decadal periods relating to the shift in surface temperature regimes (2000-2009 & 2010-2019). Growth in fishes is commonly modeled using the "Von-Bertalanffy" growth function (VBGF) to capture how a fishes size (length or weight) changes with age. For each species the growth was modeled by fitting the VGBF to the length distribution-at-age.$$L_t~=~L_{\infty}(1-e^{-K(t-t_0)})$$ Where $L_t$ is the length in cm of a species at age $t$. $L_\infty$ is the asymptotic maximum length, $K$ is \_\_\_\_\_, and $t_0$ is x intercept, or the hypothetical age at 0 cm length.```{r}# Set the function for nls() to solvevb <-vbFuns(param="Typical") ``````{r vonbert test example}#| eval: false#| echo: false#### What a single model fitting looks like ##### Datatest_dat <-filter(vonbert_species_agedat, comname =="haddock")# Get starting points from all datatest_starts <-vbStarts(length_cm ~ age, data = test_dat)# Fit the von berttest_vbert_fit <-nls( length_cm ~vb(age, Linf, K, t0),data = test_dat,start = test_starts)# Access parameters with coef()test_vbert_coef <-as.data.frame(t(coef(test_vbert_fit)))# Get Summary Info (For convergence and std. error)t_sum <-summary(test_vbert_fit)# Parameters w/ Standard Errort_params <-rownames_to_column(as.data.frame(t_sum$parameters), var ="parameter")# Convergence Infot_sum$convInfo$isConv``````{r vonbert model fitter}# Function to Pull Vonbert Coefficients, # In development: boot strapped intervalsvonbert_coef <-function(length_age_dat, start_points, min_obs =20){#### Von Bert Fitting Using: {FSA} ##### Don't run on fewer than a minimum number of observations num_aged <-nrow(length_age_dat)# If number aged less than the minimum observations:# fill output with NA's, but preserve structure to build tableif(num_aged < min_obs){# dataframe structure for coefficients na_coef_df <-data.frame("Linf"=NA,"K"=NA,"t0"=NA)# dataframe structure for parameters na_param_df <-data.frame("parameter"=NA,"estimate"=NA,"std_error"=NA,"t_value"=NA,"pr_t"=NA,"converged"=NA)# What to return when it fails/aborts:return(list("model"=NA, "coef"= na_df, "params"= na_param_df)) } # Close conditional for sample size#### Fitting Model When Sample Size is Sufficient ##### Use nls to estimate VBGF parameters using starting points vbert_fit <-nls(length_cm ~vb(age, Linf, K, t0),data = length_age_dat,start = start_points)# Access parameters with coef() vbert_coef <-as.data.frame(t(coef(vbert_fit)))# Get Summary Info (For convergence and std. error) vbert_summ <-summary(vbert_fit)# Convergence Info vbert_conv <- vbert_summ$convInfo$isConv# Parameters w/ Standard Error vbert_params <-as.data.frame(vbert_summ$parameters) %>%rownames_to_column(var ="parameter") %>% janitor::clean_names()# add convergence information in vbert_params$converged <- vbert_conv#### Work in Progress / Debugging ##### This section fails, can't find namespace/object "length_age_dat"# # Use bootstrapping for confidence intervals:# vbert_boot <- Boot(vbert_fit) # Be patient! Be aware of some non-convergence# vbert_conf <- confint(vbert_boot)# output list out_l <-list("model"= vbert_fit,"coef"= vbert_coef,"params"= vbert_params#,#"conf" = vbert_conf )# Return the listreturn(out_l)}``````{r vonbert group wrapper function}# Second Function for running that for a single Species, split by a factor column# Function does these things:# 1. grab data for one species# 2. Get starting points for VB using all data on that species# 3. Fit Vonbert to both Regimesprocess_group_vbert <-function(spec_name, split_col, min_obs){# Get the data using the species spec_dat <- dplyr::filter(vonbert_species_agedat, comname == spec_name)# Get starting points from all data test_starts <-vbStarts(length_cm ~ age, data = spec_dat)# Get coefficients for groups spec_dat %>%split(.[split_col]) %>%map(.f = vonbert_coef,start_points = test_starts, min_obs = min_obs) %>%# map_dfr(pluck, "coef", .id = split_col) %>% # Pull out just the coefficientsmap_dfr(pluck, "params", .id = split_col) %>%# Pull out the parameters and stderrmutate(comname = spec_name)}``````{r vonbert coefficient processing}# Running Each Species * Regimespecies_coef <- vonbert_species %>%map_dfr(function(x){process_group_vbert(spec_name = x, split_col ="regime", min_obs =30) })```### Growth Impacts on ProductivityYield-per-recruit analysis similar to Baudron et al. 2014?# Results## Temperature Regime Shift```{r}#| label: Temperature regime shift test# Pull the whole area's tempssurvey_temps <-filter(temp_regimes, survey_area =="Full Region")# Pull the regime averages outregime_avgs <- survey_temps %>%distinct(regime, avg_temp_c, avg_anom_c, temp_label) # Get the regime averages separately for in text valuesregime_1 <- regime_avgs %>%filter(regime =="1982-2009 Average")regime_2 <- regime_avgs %>%filter(regime =="2010-2021 Average")# Functions developed in 10_change_points.Rlibrary(bcp)# run bcp and flag changepoints using thresholdflag_bcp <-function(response_var, predictor_var =NULL, id_var, burnin =10000,mcmc =10000, return.mcmc = T, threshold_prob =0.3){# Use bcp to find support for change points bcp_x <-bcp(y = response_var, x = predictor_var, burnin = burnin, mcmc = mcmc, return.mcmc = return.mcmc)sink("/dev/null")# Convert summary to a dataframe bcp_sum <-as.data.frame(summary(bcp_x))sink()# Label the id variable bcp_sum$id <- id_var#Only pull records with support above desired threshold flag_dat <- bcp_sum[which(bcp_x$posterior.prob > threshold_prob), ] flag_dat <-as.data.frame(flag_dat)return(flag_dat)}# Set seedset.seed(3905)# Perform BCP regions temp breakpointsanom_bcp <-flag_bcp(response_var = survey_temps$sst_anom, id_var = survey_temps$yr, threshold_prob =0.3)# in-text valuesprob_text <-paste(round(anom_bcp$Probability *100, 2), collapse =", ")```SST anomaly time series of the region indicate a jump in annual temperature anomalies indicative of a regime shift centered around 2010. A bayesian change point analysis showed highest support for a regime shift among `r nrow(anom_bcp)` the years: 2009-2011. Change point probabilities for these years were `r prob_text`% respectively.```{r}#| label: bcp plotting function# Function to plot bcp breaks, highlighting their supportplot_bcp_breaks <-function(bcp_data, bcp_breaks, region_label, y_label ="sizeSpectra Exponent (b)",y_col ="b"){# Set label elevation label_y <-mean(as.data.frame(bcp_data)[, y_col], na.rm = T) -1# plot with breakpointggplot(bcp_data, aes_string("yr", y_col)) +geom_line() +geom_point() +geom_vline(data = bcp_breaks, aes(xintercept = id), linetype =2, color ="gray50") + ggrepel::geom_label_repel(data = bcp_breaks, aes(x = id, y = label_y, label =str_c(round(Probability, 2), "%"))) +labs(x ="Year", y = y_label, title =str_c(region_label, " - Bayesian Change Points"))}# Plot the bcp breaks for the Sea surface temperature anomaliesplot_bcp_breaks(bcp_data = survey_temps, bcp_breaks = anom_bcp, y_label ="Temperature Anomalies", region_label ="Northeastern US Continental Shelf", y_col ="sst_anom")```### Warming Trends```{r}#| label: temperature warming trends# warming trendslist("1982-2019"=filter(temp_regimes, between(yr, 1982, 2009)),"2010-2019"=filter(temp_regimes, between(yr, 2010, 2019))) %>%map_dfr(function(.x){ mod_x <-lm(sst ~ yr, data = .x) rate <-coef(mod_x)[["yr"]] *10 rate_df <-tibble("Decadal Rate"=round(rate, 2))return(rate_df)}, .id ="Period")```Temperatures prior to the 2010 regime shift rose at a rate of 0.7 $°C/Decade$. Regional temperatures then peaked in 2012 and have remained high while decreasing slowly at a rate of 0.02 $°C/Decade$. Average temperatures following the regime shift (2011-2021), were on average 1.09°C warmer than the than the average over the previous 28 years (1982-2009). Temperatures within each regime were relatively stable, indicative of a jump in mean temperatures rather than a rate change between them. While temperatures in the current regime are falling, the rate is such that```{r temp regimes}#| label: temperature regime figure# Plot how each one exhibits a break in the temperature across the regiontemp_regimes %>%filter(survey_area =="Full Region") %>%ggplot( aes(yr, sst_anom)) +geom_col(aes(fill = regime), size =0.75, alpha =0.4) +scale_fill_gmri() +# Using geomtextpath for labelgeom_textpath(aes(x = yr, y = avg_anom_c, label = temp_label, linecolor = regime,colour = regime), size =4, vjust =-1, fontface ="bold", show.legend =FALSE) +# Line segments for the regime averages:geom_segment(aes(x = x, xend = xend,y = avg_anom_c, yend = avg_anom_c,color = regime),#key_glyph = "rect",size = .8,linetype =1) +scale_x_continuous(expand =expansion(add =c(0.25,0.25))) +scale_y_continuous(labels =number_format(suffix =" \u00b0C")) +scale_color_gmri() +guides(fill ="none") +theme(legend.title =element_blank(),legend.position ="bottom",legend.background =element_rect(fill ="transparent"),legend.key =element_rect(fill ="transparent", color ="transparent")) +labs(title ="Northwest-Atlantic Temperature Regime Shift",x ="", y ="Surface Temperature Anomaly",caption ="Anomalies calculated using 1982-2011 reference period.") +guides(label ="none",color =guide_legend(keyheight =unit(0.5, "cm"), ))```## Distribution Shifts```{r}#| fig.height: 8#| fig-width: 8# Get the values for each of them in a standardized format# Standardize to global benchmarks# Latitudelat_yrly <-biomass_weighted_var(station_totals, decdeg_beglat, est_year) %>%left_join(select(global_benchmarks, c("comname", "avg_lat", "lat_sd")), by ="comname") %>%mutate(biom_z = (biom_mu - avg_lat) / lat_sd) # Longitudelon_yrly <-biomass_weighted_var(station_totals, decdeg_beglon, est_year) %>%left_join(select(global_benchmarks, c("comname", "avg_lon", "lon_sd")), by ="comname") %>%mutate(biom_z = (biom_mu - avg_lon) / lon_sd) # Bottom Temperaturet_yrly <-biomass_weighted_var(station_totals, bottemp, est_year) %>%left_join(select(global_benchmarks, c("comname", "avg_temp", "temp_sd")), by ="comname") %>%mutate(biom_z = (biom_mu - avg_temp) / temp_sd) # Depthd_yrly <-biomass_weighted_var(station_totals, avgdepth, est_year) %>%left_join(select(global_benchmarks, c("comname", "avg_depth", "depth_sd")), by ="comname") %>%mutate(biom_z = (biom_mu - avg_depth) / depth_sd) # Put them together in one table for heatmapstandardized_changes <-list("Latitude"= lat_yrly,"Longitude"= lon_yrly,"Bottom Temperature"= t_yrly,"Depth"= d_yrly) %>%map_dfr(function(x){select(x, comname, est_year, var, biom_z)}, .id ="var_neat")# Organize them by their global average latitudeglobal_benchmarks <- global_benchmarks %>%mutate(comname =fct_reorder(.f = comname, .x = avg_lat, .fun = median))# Relevel the standardized changesstandardized_changes <- standardized_changes %>%mutate(comname =factor(comname, levels =levels(global_benchmarks$comname)))# Heatmapstandardized_changes %>%ggplot(aes(est_year, comname, fill = biom_z)) +geom_tile() +scale_x_continuous(expand =expansion(add =c(0, 0))) +facet_wrap(~var_neat, ncol =2, scales ="fixed") +scale_fill_distiller(palette ="RdBu",oob = scales::oob_squish,limits =c(-1, 1), breaks =c(-1, -.5, 0, .5, 1),labels =c("<-1", "-.5", "0", ".5", ">1")) +labs(x ="Year", y ="",title ="Center of Biomass Movements Through Time",fill ="Center of Biomass Metric Change (z)",caption ="Species ordered based on center of biomass latitude across all years.") +theme(legend.position ="bottom",axis.text.y =element_text(size =7)) +guides(fill =guide_colorbar(title.position ="top", title.hjust =0.5, barwidth =unit(6, "cm")))```## Center of Biomass Trends::: panel-tabset### Latitude```{r}# Standardized Plotslat_yrly %>%ggplot(aes(est_year, biom_z, group = comname)) +# geom_path(aes(group = comname), linetype = 1, color = "gray50", size = 0.5, alpha = 0.2) +geom_smooth(method ="lm", formula = y ~ x, alpha =0.8, se =FALSE, color ="gray50") +geom_jitter(width =0.1, size =2, alpha =0.25) +geom_hline(aes(yintercept =0, color ="Overall Average"), size =1, linetype =1) +scale_color_gmri(reverse = T) +labs(title ="Shift in Center of Biomass Latitudes",subtitle ="Yearly averages weighted by biomass, scaled against all-year average.",y ="Center of Biomass Latitude (z)", x ="", color ="Time Period", fill ="Time Period") +theme(legend.position ="bottom")```### Longitude```{r}# Longitudelon_yrly %>%ggplot(aes(est_year, biom_z, group = comname)) +# geom_path(aes(group = comname), linetype = 1, color = "gray50", size = 0.5, alpha = 0.2) +geom_smooth(method ="lm", formula = y ~ x, alpha =0.8, se =FALSE, color ="gray50") +geom_jitter(width =0.1, size =2, alpha =0.25) +geom_hline(aes(yintercept =0, color ="Overall Average"), size =1, linetype =1) +scale_color_gmri(reverse = T) +#coord_flip() +labs(title ="Shift in Center of Biomass Longtitudes",subtitle ="Yearly averages weighted by biomass, scaled against all-year average.",y ="Center of Biomass Longitude (z)", x ="", color ="Time Period", fill ="Time Period") +theme(strip.text.y =element_text(size =7, angle =0, hjust =0),legend.position ="bottom")```### Mean Temperature```{r}# Temperaturet_yrly %>%ggplot(aes(est_year, biom_z, group = comname)) +# geom_path(aes(group = comname), linetype = 1, color = "gray50", size = 0.5, alpha = 0.2) +geom_smooth(method ="lm", formula = y ~ x, alpha =0.8, se =FALSE, color ="gray50") +geom_jitter(width =0.1, size =2, alpha =0.25) +geom_hline(aes(yintercept =0, color ="Overall Average"), size =1, linetype =1) +scale_color_gmri(reverse = T) +#coord_flip() +labs(title ="Shift in Center of Biomass Temperature",subtitle ="Yearly averages weighted by biomass, scaled against all-year average.",y ="Center of Biomass Bottom Temperature (z)", x ="", color ="Time Period", fill ="Time Period") +theme(legend.position ="bottom")```### Mean Depth```{r}# Depthd_yrly %>%ggplot(aes(est_year, biom_z, group = comname)) +# geom_path(aes(group = comname), linetype = 1, color = "gray50", size = 0.5, alpha = 0.2) +geom_smooth(method ="lm", formula = y ~ x, alpha =0.8, se =FALSE, color ="gray50") +geom_jitter(width =0.1, size =2, alpha =0.25) +geom_hline(aes(yintercept =0, color ="Overall Average"), size =1, linetype =1) +scale_color_gmri(reverse = T) +#coord_flip() +labs(title ="Shift in Center of Biomass Depth",subtitle ="Yearly averages weighted by biomass, scaled against all-year average.",y ="Center of Biomass Depth (z)", x ="", color ="Time Period", fill ="Time Period") +theme(strip.text.y =element_text(size =7, angle =0, hjust =0),legend.position ="bottom")```:::## Changes to Growth```{r regime comparison data}# Take the first part from the previous plot:# takes each species# for age 0 - max_age, solve what the size would be based on vb fitregime_curve_prep <-function(spec_name){# The raw size/age data: t_dat <- vonbert_species_agedat %>%filter(comname == spec_name) # The group coefficients, augmented with x values across a range age_max <-max(t_dat$age) x_range <-seq(from =0, to = age_max, by = .1) t_coef <- species_coef %>%distinct(comname, regime, parameter, estimate) %>%pivot_wider(values_from ="estimate", names_from ="parameter") %>%filter(comname == spec_name) %>%right_join(# Make a dataframe that covers the size range for the speciesy =data.frame(regime =rep(unique(t_dat$regime), length(x_range)),age =rep(x_range, length(unique(t_dat$regime)))),by ="regime") %>%mutate(# Calculate the length at agepred_length = Linf * (1-exp(-1* K * (age - t0))) )# Pivot wider?# Take the size at age from each regime and get the difference t_coef <- t_coef %>%select(-c(Linf, K, t0)) %>%pivot_wider(names_from ="regime", values_from ="pred_length") %>%mutate(growth_shift =`Warm Regime`-`Early Regime`,growth_frac = (growth_shift /`Early Regime`),size_change =ifelse(growth_shift >0, "Increase", "Decrease") )return(t_coef)}# Do it for all the speciesregime_comparison <- vonbert_species %>%map_dfr(.f = regime_curve_prep, .id ="comname")```Changes in the size-at-age growth relationship varied across species, with the majority (14/18) of the species exhibiting smaller L-infinity during the warmer temperature regime.```{r}#| eval = TRUE,#| fig.height = 2# Try geom_waffle to just count each species as an individual and show that they are growing to smaller sizes in new regime:library(ggwaffle)age_set <-4ggplot(data =waffle_iron(data =filter(regime_comparison, age == age_set), aes_d(group = size_change), rows =2), aes(x, y, fill =fct_rev(group))) +geom_waffle() +theme_waffle() +scale_fill_gmri() +theme(axis.text =element_blank(),axis.text.x =element_blank(),axis.text.y =element_blank(),legend.title =element_text(face ="bold", size =12),legend.position ="bottom",panel.grid.major.y =element_blank()) +labs(x ="", y ="", fill =str_c("Warm Regime Growth Impacts on ", age_set," Body Length:"))```The species that showed larger asymptotic lengths include the three hake species: **red hake, silver hake, & white hake** as well as **Atlantic herring.** These species are all among the majority group of species that is found primarily in the Northern part of the survey area, in and around the Gulf of Maine.```{r, fig.height = 6}#| eval = TRUE# Reshape things so you can show confidence intervals hereplot_dat <- species_coef %>%filter(parameter !="t0") %>%mutate(llim = estimate -2* std_error,rlim = estimate +2* std_error,residence =ifelse( comname %in%c("atlantic mackerel", "scup", "black sea bass", "butterfish"),"Southern\nResident", "Gulf of Maine\nResident"),regime =fct_rev(regime),comname =fct_rev(comname))# Plot itggplot(plot_dat, aes(x = estimate, y = comname, group = regime, color = regime)) +# need to nudge groups separately# Early Regimegeom_segment(data =filter(plot_dat, regime =="Early Regime"), aes(x = llim, xend = rlim, yend = comname), position =position_nudge(x =0, y =0.25)) +geom_segment(data=filter(plot_dat, regime =="Warm Regime"), aes(x = llim, xend = rlim, yend = comname), position =position_nudge(x =0, y =-0.25)) +geom_point(aes(color = regime), position =position_dodge(width =1), size =2) +#facet_grid(. ~ parameter, scales = "free", space = "free_y") +facet_grid(residence ~ parameter, scales ="free", space ="free_y") +scale_color_gmri(reverse = F) +labs(x ="Coefficient Estimate", y ="", fill ="Temperature Regime:") +theme(legend.position ="bottom", panel.background =element_rect(color ="gray80", fill ="transparent", size =1),panel.spacing =unit(2, "lines"),panel.border =element_rect(color ="gray80", fill ="transparent", size =1))``````{r, fig.height = 6}#| eval = FALSE# Plotting Everythingspecies_coef %>%distinct(comname, regime, parameter, estimate) %>%pivot_wider(values_from =c("estimate"), names_from ="parameter") %>%pivot_longer(names_to ="vbcoef", values_to ="val", cols =c("K", "Linf", "t0")) %>%filter(vbcoef !="t0") %>%mutate(residence =ifelse(comname %in%c("atlantic mackerel", "scup", "black sea bass", "butterfish"), "Southern\nResident", "Gulf of Maine\nResident")) %>%ggplot(aes(x = val, y =fct_rev(comname), fill =fct_rev(regime))) +geom_col(position ="dodge") +# facet_grid(. ~ vbcoef, scales = "free", space = "free_y") +facet_grid(residence ~ vbcoef, scales ="free", space ="free_y") +scale_fill_gmri(reverse = T) +labs(x ="Coefficient Estimate", y ="", fill ="Temperature Regime:") +theme(legend.position ="bottom", panel.background =element_rect(color ="gray80", fill ="transparent", size =1),panel.spacing =unit(2, "lines"),panel.border =element_rect(color ="gray80", fill ="transparent", size =1))```## Changes in $L_\infty$```{r}species_coef %>%distinct(comname, regime, parameter, estimate) %>%pivot_wider(values_from ="estimate", names_from ="parameter") %>%mutate(regime_2 =ifelse(regime =="Early Regime", "early", "warm")) %>%select(comname, regime_2, Linf) %>%pivot_wider(values_from = Linf, names_from = regime_2) %>%mutate(linf_diff = warm - early,diff_perc = (linf_diff / early) *100) %>%mutate(across(where(is.numeric), round, 2)) %>%arrange(desc(diff_perc)) %>%select(c(1, 4, 5)) %>%setNames(c("Common Name", #"Early Regime Linf", "Warm Regime Linf", "Change in Linf (cm)", "Percent Change")) %>%gt()```## VB Fits to Data:```{r}# Make the plot for each speciesvb_fits_plot =function(spec_name){# The raw size/age data: t_dat <- vonbert_species_agedat %>%filter(comname == spec_name) # The group coefficients, augmented with x values across a range age_max <-max(t_dat$age) x_range <-seq(from =0, to = age_max, by = .1) t_coef <- species_coef %>%distinct(comname, regime, parameter, estimate) %>%pivot_wider(values_from ="estimate", names_from ="parameter") %>%filter(comname == spec_name) %>%right_join(y =data.frame(regime =rep(unique(t_dat$regime), length(x_range)),age =rep(x_range, length(unique(t_dat$regime)))),by ="regime") %>%mutate(# Calculate the length at agepred_length = Linf * (1-exp(-1* K * (age - t0))) )# Make plot?ggplot(t_dat, aes(x = age, y = length_cm)) +geom_point(alpha =0.25) +#facet_wrap(~regime, ncol = 1) +geom_line(data = t_coef,aes(x = age, y = pred_length, group = regime, color = regime), size =1) +scale_x_continuous(breaks =seq(0, age_max, by =2),limits =c(0, age_max),expand =expansion(mult =c(0.1, .2))) +scale_color_gmri() +geom_dl(data = t_coef, aes(x = age, y = pred_length, label = regime, color = regime), method=list("last.qp", cex =0.65, hjust =-.15),size =12) +labs(x ="Age",y ="Length (cm)",title = spec_name,color ="Temperature Regime:",subtitle ="Length at Age Regimes") }```::: panelset```{r, results="asis", fig.height = 6}# Loop through the key species:walk(key_species, plot_panelset, plot_fun = vb_fits_plot)```:::# Discussion------------------------------------------------------------------------# Bonus: Data Explorations## Age Data AvailabilityThe availability of age data varied by species and through time, with older ages becoming less frequently sampled in recent years across many species.```{r ridgeplot fun}# Plotting age distribution as ggridgesplot_age_ridgeplot <-function(species){ vonbert_species_agedat %>%filter(comname == species, age <=15) %>%mutate(yr =factor(est_year),yr =fct_rev(yr)) %>%ggplot(aes(x = age, y = yr)) +geom_density_ridges(aes(fill = regime), color ="white", show.legend =FALSE) +scale_fill_gmri() +facet_wrap(~comname) +scale_x_continuous(limits =c(0,NA), expand =expansion(add =c(0, 0))) +labs(x ="Age", y ="Year") +guides(size =guide_legend(title.position ="top", title.hjust =0.5))}```::: panelset```{r ridge panels, results="asis", fig.height=8}# Loop through the key species:walk(key_species, plot_panelset, plot_fun = plot_age_ridgeplot)```:::Strong age classes were particularly prominent in the following species, showing lasting numbers in subsequent years.| Species | Strong Year Classes ||------------------|---------------------|| Atlantic Cod | 2003, 2013 || Haddock | 2003, 2010, 2013 || Atlantic Herring | 2008, 2011 || Pollock | 2000, 2012 |## Growth Increment ChangesAnnual Growth increments for each species were computed to track the change in growth of specific age-groups through time.```{r}# Age-Size Increments Differencesvb_incs <- vonbert_species_agedat %>%group_by(`common name`, est_year, age, regime) %>%summarise(avg_len =mean(length_cm, na.rm = T),avg_wt =mean(indwt, na.rm = T),.groups ="drop") %>%arrange(`common name`, age) %>%group_split(`common name`, est_year) %>%map_dfr(function(x){ x %>%mutate(# What age is being subtracted fromlead_age =lead(age, n =1),# String built from the value to confirmdiff_name =str_c(age, " - ", lead_age),# Length/Weight at age + 1len_age_p1 =lead(avg_len, n =1),wt_age_p1 =lead(avg_wt, n =1),# Length/Weight Growth Incrementslen_inc = len_age_p1 - avg_len,wt_inc = wt_age_p1 - avg_wt,# Percent change in bodymass for going from age to age+1len_perc_change = (len_inc / avg_len) *100,wt_perc_change = (wt_inc / avg_wt) *100) %>%filter((lead_age - age ==1)) })```Age specific growth (annual growth increments) help detail where in a species life cycle it may be experiencing favorable or less-favorable conditions for growth. These can be a result of differences in habitat use or prey resources at different ages, and/or a change in exposure to temperature or other stressors that inflict a metabolic cost.```{r}# Put together a timeseries of:# how big were the age 0's,# how big were the age 1's# Plotting Age Incrementsplot_age_increments <-function( spec_choice, var_choice ="length", gam_knots =15, facet_prefix ="Age Increment: ",facet_suffix =""){# Code to switch variables and their labels: var_sym <-switch(var_choice,"length"=sym("len_inc"),"weight"=sym("wt_inc")) var_label <-switch(var_choice,"length"="Annual Growth (cm)","weight"="Annual Growth (kg)")# Add pre/suffix to facet labels appender <-function(string, prefix = facet_prefix, suffix = facet_suffix) {paste0(prefix, string)}# Make the plot vb_incs %>%filter(`common name`== spec_choice, lead_age <=5) %>%drop_na({{var_sym}}) %>%ggplot(aes(est_year, {{var_sym}})) +# geom_col(fill = gmri_cols("gmri blue")) +geom_point(color =gmri_cols("gmri blue")) +geom_line(color =gmri_cols("gmri blue"), group =1, alpha =0.3, size =0.8) +geom_smooth(method ="gam",formula = y ~s(x, bs ="cs", k = gam_knots),se =FALSE,size =1,color =gmri_cols("orange")) +facet_wrap(~diff_name, ncol =1, scales ="free",labeller =as_labeller(appender)) +theme(legend.position ="bottom") +labs(x ="",y = var_label,color ="Age Increment",title =str_to_title(spec_choice),subtitle =str_c("Changes in annual ", var_choice, " growth increments."))}```::: panelset```{r, results="asis", fig.height = 8}# Loop through the key species:walk(key_species, plot_panelset, plot_fun = plot_age_increments, var_choice ="weight", gam_knots =5)```:::## Length at Age Changes```{r regime size plots}# Distribution of Length/Weight Across Regimes plot_regime_sizes <-function(spec, var ="length", age_lim =5, facet_prefix ="Age - ", facet_suffix =""){# switches for length or width var_col <-switch (var,"length"=sym("length_cm"),"weight"=sym("indwt")) var_lab <-switch (var,"length"="Individual Length (cm)","weight"="Individual Weight (kg)")# Add pre/suffix to facet labels appender <-function(string, prefix = facet_prefix, suffix = facet_suffix) {paste0(prefix, string)}# Select Ages, Prep data size_changes <- regime_dat %>%filter( comname == spec, age <=6) %>%## Calculate group averagesgroup_by(age, regime) %>%mutate(median_size =median(!!{{var_col}})) %>%ungroup() # Make median dataframe median_df <- size_changes %>%distinct(age, regime, median_size) %>%pivot_wider(names_from ="regime", values_from ="median_size")# Make Plot size_changes %>%ggplot(aes(x =!!{{var_col}})) +## distributionstat_halfeye(aes(color = regime, fill =after_scale(color)),slab_alpha = .5, .width =0, trim = T, shape =21, point_colour ="gray25",stroke =1.6,scale = .86) +# median linegeom_segment(data = median_df,aes(x =`Early Regime`,xend =`Warm Regime`,y =0,yend =0),size = .8,color ="grey25", # , position = position_nudge(y = -.15) ) +## median points - not workingstat_halfeye(aes(color = regime), .width =c(0), slab_fill =NA) +facet_wrap(~age, ncol =1, labeller =as_labeller(appender), scales ="free") +scale_color_gmri() +scale_fill_gmri() +scale_y_continuous(labels =percent_format()) +theme(legend.position ="bottom") +labs(x = var_lab,y ="Proportion",title = spec,fill ="Temperature Regime",color ="Temperature Regime",shape ="Median" )}# Tester: # plot_regime_sizes(key_species[[1]], var = "length")```::: panelset```{r, results="asis", fig.height = 8}# Loop through the key species:walk(key_species, plot_panelset, plot_fun = plot_regime_sizes, var ="length")```:::## Weight at Age Changes::: panelset```{r, results="asis", fig.height = 8}# Loop through the key species:walk(key_species, plot_panelset, plot_fun = plot_regime_sizes, var ="weight")```:::## Shifts in Growth```{r regime impact plot}#| fig.height = 8,#| fig.width = 8# Plot their growth curvesregime_comparison %>%#filter(comname %in% key_species) %>% ggplot(aes(x = age, y = growth_frac)) +geom_hline(yintercept =0, size =0.5, lty =2, color ="black") +geom_line(aes(group = comname), show.legend =FALSE, size =1, color ="orange")+scale_y_continuous(labels =percent_format()) +facet_wrap(~comname, scales ="free_x", ncol =4) +labs(x ="Age", y ="Shift in Length Moving to Warm Regime", title ="Thermal Regime Impacts on Growth",subtitle ="Changes in Size at Age (Warm Regime - Historic Regime)") ``````{r}#| eval = FALSE# Kathy requested the curve data.# Going to save both the coefficients and the fitted data for all specieslibrary(here)write_csv(regime_comparison, here("data/growth_results/growth_regime_fits.csv"))write_csv(species_coef, here("data/growth_results/growth_regime_coefs.csv"))``````{r}# Try horizonslibrary(ggHoriPlot)# Filter outliers beyond 5th and 95th quantilecutpoints <- regime_comparison %>%mutate(outlier =between( growth_frac,quantile(growth_frac, 0.05, na.rm = T)-1.5*IQR(growth_frac, na.rm = T),quantile(growth_frac, 0.95, na.rm = T)+1.5*IQR(growth_frac, na.rm = T))) %>%filter(outlier)# Set originori <-0# Manually pick cut pointssca <-seq(-1, 1, length.out =9)sca_bins <-str_c(sca[1:(length(sca)-1)], " to ", sca[2:length(sca)], "%")sca_labels <-rev(sca_bins)theme_set(theme_bw() +theme(# Titlesplot.title =element_text(hjust =0, face ="bold", size =14),plot.subtitle =element_text(size =9),plot.caption =element_text(size =8, margin =margin(t =20), color ="gray40"),plot.margin =unit(c(1, 1, 2, 1), "lines"),legend.title =element_text(size =9),legend.text =element_text(size =9),# Axesaxis.line.y =element_line(color ="black"),axis.ticks.y =element_line(), axis.line.x =element_line(color ="black"),axis.ticks.x =element_line(), axis.text =element_text(size =11),axis.title =element_text(size =12),rect =element_rect(fill ="transparent", color ="black"),# Facetsstrip.text =element_text(color ="white", face ="bold",size =11),strip.background =element_rect(color ="#00736D", fill ="#00736D", size =1, linetype="solid")))# Plot it allggplot(regime_comparison, aes(x = age, y = growth_frac)) +geom_horizon(aes(age, y = growth_frac,fill = ..Cutpoints..), origin = ori, horizonscale = sca) +scale_fill_hcl(palette ='RdBu', reverse = T, labels = sca_labels) +facet_grid(comname~.) +theme_bw() +scale_x_continuous(limits =c(0,10), expand =expansion(add =c(0, 0))) +theme(# Titlesplot.title =element_text(hjust =0, face ="bold", size =14),plot.subtitle =element_text(size =9),plot.caption =element_text(size =8, margin =margin(t =20), color ="gray40"),plot.margin =unit(c(1, 1, 2, 1), "lines"),legend.title =element_text(size =9),legend.text =element_text(size =9),# Axesaxis.line.y =element_line(color ="black"),axis.line.x =element_line(color ="black"),axis.ticks.x =element_line(), axis.text =element_text(size =11),axis.title =element_text(size =12),rect =element_rect(fill ="transparent", color ="black"),# Facetsstrip.text =element_text(color ="white", face ="bold",size =11),strip.background =element_rect(color ="#00736D", fill ="#00736D", size =1, linetype="solid"),panel.grid =element_blank(),panel.spacing.y=unit(0.05, "lines"),strip.text.y =element_text(size =7, angle =0, hjust =0, face ="bold"),axis.text.y =element_blank(),axis.title.y =element_blank(),axis.ticks.y =element_blank(),panel.border =element_blank(),legend.position ="left") +labs(x ='Age', fill ="Change in Body Size at Age") +guides(fill =guide_legend(ncol =1))````r insert_gmri_footer()`
